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We present a generalization of the Green-Kubo expressions for thermal transport coefficients \i in 
complex fluids of the generic form, /j, = fi x + dtV~ x (J e exp(t£) J)o, i.e. a sum of an instantaneous 
transport coefficient /loo, and a time integral over a time correlation function in a state of thermal 
equilibrium between a current J and its conjugate current J e . The streaming operator exp(f£) 
generates the trajectory of a dynamical variable J(t) — exp(t£)J when used inside the thermal 
average (• • ■ )o- These formulas are valid for conservative, impulsive (hard spheres), stochastic and 
dissipative forces (Langevin fluids), provided the system approaches a thermal equilibrium state. 
In general /i^ and J e 7^ J, except for the case of conservative forces, where the equality signs 
apply. The most important application in the present paper is the hard sphere fluid. 

PACS numbers: 05.40.-a, 05.20.Dd 



As the interest in fluids during the last decades has been shifting from standard fluids with smooth conservative 
interactions to more complex fluids with ditto interactions, a new analysis of linear response theory for such systems 
seems timely. The present paper provides this analysis for a wide class of complex fluid-type systems. 

The Green-Kubo formulas for thermal transport coefficients in simple classical fluids with conservative interactions 
are widely used, and generally accepted 0,0,1110] as exact expressions for general densities, as long as the deviations 
from equilibrium and the gradients are small, and the transport coefficients exist. These expressions are given in 
terms of equilibrium time correlation functions between ./V-particle currents, i.e. 



where /1 denotes a typical transport coefficient, (. . . }o is an average over a thermal equilibrium ensemble at temperature 
T = I/UbP, with a weight po ~ exp[— 0H], and limy^oo denotes the thermodynamic limit. In the sequel this limit is 
understood, but not explicitly written. The time evolution of a dynamical variable, J(t) = e tc J(0), can be described 
by a streaming operator e t£ , and is invariant under the time reversal transformation. The infinitesimal generator, 
C, referred to as Liouville operator, changes sign under the time reversal transformation. The total microscopic 
flux J = linik^o jk is related to the Fourier mode Ok of a conserved density through the local conservation law, 
(9 t ak = £a,k = — i\t • jk, and contains in general contributions from kinetic transport, and from collisional transfer, 
i.e. transport through inter-particle forces. The simple form II. 1|) only applies to systems with time reversal invariant 
dynamics, as we shall see in this paper. 

A summary of the results for systems with non-conservative forces has been given in Ref. 0] , and the goal of this 
article is to provide a general analysis with detailed applications to hard sphere fluids. The applications to stochastic 
systems, as described by Langevin and Fokker-Planck equations will be given elsewhere 6]. Green-Kubo formulas 
still exist, but their generic form is different from l|I.ljl . The systems of interest are elastic hard sphere fluids with 
impulsive forces Q, and systems with Brownian dynamics lJ3i considered as a mixture of hard spheres where the 
mass ratio is taken to infinity; lattice gas cellular automata [!j and multi-particle collision dynamics (MPCD) models, 
where the discrete time dynamics involves stochastic variables 0. Finally there is the large class of complex fluids, 
where TV-particle fluid or magnetic systems are described by mesoscopic Langevin equations contai ning dissipative 
and stochastic for ces, or by corresponding Fokker-Planck equations for the probability distribution |lll H^ITii. HH 
13 13 IH 13 13 Efl El Hi 13] In fact, also the models for critical and unstable fluids and magnets Hj belong 
to the Langevin models covered by the present theory. We also note that the Langevin and Fokker-Planck equations 
are not time reversal invariant. The theory to be developed here may be applied as soon as the analytic form of the 
generator C and of the thermal distribution po are given, i.e. as soon as the time evolution can be expressed either in 
terms of hard sphere dynamics or in terms of Langevin and Fokker-Planck equations. 

The systems listed above not only refer to simple fluids, but cover a large range of complex fluids, also outside the 
collection of critical and unstable systems. For instance, Ref. [3] considers diffusion coefficients of colloidal particles, 



I. INTRODUCTION 




(LI) 



whereas Ref. (23 contains an extensive list of applications of the so-called dissipative particle dynamics (DPD) fluid 
to colloidal particles, emulsions and polymers |lllll2j|. These are just a few examples, but there are many more non- 
equilibrium fluid- type models, that take fluctuations into account, satisfy detailed balance criteria, and reach thermal 
equilibrium, for which Green-Kubo formulas have been or can be derived. 

Granular fluids are outside the scope of this study, although Green-Kubo formulas have been considered in the 
literature [25l l2fi| . The reason is that such systems do not approach a stable state of thermal equilibrium, which is 
a basic requirement in the present derivation. In the above articles the stable equilibrium state is replaced by the 
unstable homogeneous cooling state. Such extensions are outside the scope of the present paper. 

The goal of this study is to derive generalized Green-Kubo expressions for large classes of TV-particle systems, that 
do approach thermal equilibrium, and have equations of motion that possibly lack time reversal invariance. The 
interaction forces may be impulsive, dissipative, stochastic, and possibly include conservative forces as well. The 
method used is linear response theory [lj, which yield exact expressions for transport coefficients, valid for arbitrary 
densities. They have the generic form, 

/•OO 

H = Moo + V- 1 / dt(J t e tc J) . (1.2) 
Jo 

The pseudo-streaming operators exp(i£) with t > are well-defined for all points in the relevant phase space, and 
generate the proper trajectories of the dynamical variables when used inside statistical averages (• • • } . An intuitive 
generalization of (|I.lJl to non-conservative forces does not seems obvious. For the general class of Langevin systems a 
result similar to p. 2(1 has been obtained in Ref.[2j|, but the instantaneous contribution /Jqo to the dissipative transport 
coefficient has not been identified. The only known results of the complete form 111.21) with instantaneous part as well 
as time correlation part have been obtained in Ref. 8] for the friction coefficient of Brownian particles, and in Ref.Q 
for the shear viscosity of the hard sphere fluid. In the former case the second term in (11.21) contains a time correlation 
function of an inter-particle force and its conjugate force, in the latter of a stress and a conjugate stress. 

A discussion about exact expressions for transport coefficients is not complete without the Einstein-Helfand formulas 
for the Navier-Stokes transport coefficients 0, [3, Q • They are the analogs of Einstein's formula for the self diffusion 
D s , and are also expressed in equilibrium time correlation functions, e.g. 

D s = lim (l/2t) ((x(t) - x(0)) 2 ) 
n = lim ifi/2Vt) ((M(t) - A/(0)) 2 ) , (1.3) 

where n is the shear viscosity with M(t) = Yli mv iv{t) x i(t)- ^ n ^ ac ^ these forms give a more compact and more 
robust representation of transport coefficient than the Green-Kubo formulas {QJl or (|I.2|1 • This is the case because 
the former can be expressed by only specifying the momentum and energy of a particle, and do not require the explicit 
introduction of forces. Consequently, the Helfand formulas apply to systems with conservative, impulsive, dissipative 
and stochastic forces, and either type of Green-Kubo formulas (|I.I|) or i|l.2l) can be derived for the relevant cases. An 
added advantage of the Einstein-Helfand formulas is that they are evidently positive. 

So there exist two independent routes for generalizing Green-Kubo formulas to non-conservative forces. The first 
one is to start from the Einstein-Helfand formulas, and transform them into (|I.2f) . This route has been followed in 
Ref.0- The second one is to derive the generalizations from the fundamental Liouville equation, or its mesoscopic 
analog, the Fokker-Planck equation. The latter route is followed in the present paper. 

The plan of the paper is as follows: in Sec. II we obtain a fundamental commutation relation between the N— particle 
equilibrium distribution function po and the infinitesimal generator C (Liouville or Fokker-Planck operator), that is 
required in the derivation of the Green-Kubo formulas. In Sec. Ill the generic form lll.2t of these formulas is derived 
using the Zwanzig-Mori projection operator method in a form that is applicable to systems with conservative forces 
(classical fluids, lattice gas cellular automata, MPCD models), impulsive forces (hard spheres), and stochastic and 
dissipative forces (Langevin fluids). In Sec. IV we apply the new results to hard sphere fluids with some technical 
details deferred to Appendix A. Further comments and conclusions are presented in Sec.V. 

II. FUNDAMENTAL COMMUTATION RELATION 

We will analyze the similarities and differences among the systems, mentioned in the introduction, all of which 
evolve to a thermal equilibrium state described by the N— particle equilibrium distribution function pq. This enables 
us to formulate a fundamental commutation relation for the infinitesimal generator £, that is required in the derivation 
of the Green-Kubo formulas for systems whose equations of motion are not time reversal invariant, or more precisely, 



whose pseudo-streaming operators (Langevin fluids, elastic hard spheres) do not generate time reversal invariant 
dynamics. 

Suppose we have constructed for the system of interest the infinitesimal generator C and the streaming operator 
exp(tC), that generates the time evolution of dynamic variables inside averages. To study equilibrium time correlation 
functions, we need two different types of transposes of C. The first transpose, C, referred to as adjoint, is defined 
with respect to an inner product with a weight equal to unity, and a second transpose, C e , referred to as conjugate, 
is defined with a weight po, i.e., 



dTACB = J dTBCA 
(ACB) Q = J dTp ACB 

dTp BC e A = {BC e A) Q . (Ill) 



The corresponding streaming operators, A(t) = exp(tC)A(0), and A(— t) = cxp(tC e )A(0) with t > 0, generate the 
forward and backward or time reversed evolution inside averages in the sense that equilibrium time correlation functions 
are stationary, i.e. depend only on time differences, 

(A(0)5(*)>o = ( AetCB ) 
= (Be t£ A) Q = (A(-t)B(0)) Q . (11.2) 

Having obtained a conjugate generator C e , we also introduce a conjugate current, J c = limfc^oJ*k, through the 
relations £ e a£ = — ik • jj k . For conservative forces the /^-operators are first order differential operators. So, the 
adjoint C = —C. The conjugate operator C e for conservative interactions follows directly from the second line in 
(|II.lj) . because C and po commute, and the conjugate operator equals the adjoint, i.e. for time reversal invariant 
dynamics, 

C f - = -C = C. (11.3) 

Note that here the conjugate current satisfies J e = J. 

In systems, such as lattice gas cellular automata and MPCD models, the dynamics involves stochastic variables, but 
does not contain any dissipative interactions. The microscopic equations of motion inside averages are time reversal 
invariant. Here the statistical averages include averages over all realizations of the stochastic variables. The time 
reversal transform of the Liouville operator C e satisfies the relation, C e = — C, as is the case for conservative systems. 
This transform is obtained by reversing the sign of all variables odd under the time reversal transformation (e.g. 
velocities). Of course even variables remain unchanged (e.g. positions, energies). 

Next, we consider the Langevin fluids. They are described by iV-particle Langevin equations [13,123, an d reach 
thermal equilibrium. From these equations the corresponding Fokker-Planck equation for the iV-particle probability 
distribution p(T, t) is constructed, which is the fundamental starting point for our linear response theory. The latter 
is written here 29] as dtp = Cp. As we are interested in time correlation functions, it is sufficient to define streaming 
operators e tc for the time evolution of dynamical variables inside averages. This may be done through the adjoint, 

(A(t)) = J dTAe tZ p(T, 0) = J dTp(T, 0)e tc A = (e tc A). (II.4) 

In this formulation averages over the realizations of stochastic forces are already accounted for in the Fokker-Planck 
equation. So, C is the adjoint of the Fokker-Planck operator C. Because of the presence of dissipative forces, the 
Fokker-Planck equation is not time-reversal invariant. Consequently, C e ^ — C (no odd parity in the velocities). The 
Fokker-Planck operator C and its adjoint C contain first and second order derivatives, and also here the adjoint can 
be simply obtained. 

These Langevin equations have been constructed such that the detailed balance condition holds by imposing that 
the Gibbs' distribution is a steady state solution of the Fokker-Planck equation, i.e. dtpo = Cpo = 0. Then dissipative 
and random forces are related by the fluctuation-dissipation theorem. This implies that the Fokker-Planck equation 
has an //-theorem or Lyapunov functional, and that pg is the unique steady state solution. 

Next consider the construction of the conjugate operator C e , where po and C in IjlLIjl do not commute. So we 
introduce an operator C through the commutation relation 



PqC = Cp Q . 



(11.5) 



Having obtained C, we insert this relation in line 2 of 1)11. 1)) , and obtain the conjugate generator as the adjoint of C, 
i.e. 

C e = C. (11.6) 

For the Fokker-Planck equations, satisfying the detailed balance criteria, the commutation relation (|II.5(I has been 
derived in the literature 27, 28] with a general prescription for C. Only in the case of conservative interactions, where 
the equations of motion are time reversal invariant, the equalities C = C and J e — J hold. For non-conservative 
interactions one has C ^ C and J e J. 

As we will see in the next section, the replacement of 1)11.3)) for a time evolution with time reversal invariance, by 
the commutation relation 1)11. 5|l for a time evolution without, has severe effects on the explicit form of the Green- 
Kubo formulas for the Langevin fluids, as shown in 1)1. 2|) . There are two conspicuous differences with respect to 
the standard Green-Kubo formulas 1)1.1)) . The fluxes J c and J inside the equilibrium time correlation function are 
necessarily unequal, as C e ^ —C, and there appears an additional instantaneous transport coefficient p^, that solely 
depends on equilibrium properties. 

Finally we discuss fluids of elastic har d sp heres with mass m and diameter a. In the kinetic theory of hard 
sphere fluids the standard notation 0, yCJ, |3ll |33| is: C = L + , and C e = —L-. Here the pseudo-streaming operator 
exp[i£] = exp[tL + ] with t > generates forward trajectories, and exp[i£ e ] = exp[— tLJ\ with t > backward or time 
reversed trajectories. Furthermore, it has been shown [3fT | that the correct hard sphere dynamics inside an equilibrium 
average is generated by the combination poA(t) — p e tL+ A(0). Here the presence of the factor p guarantees that 
a zero weight is given to the unphysical overlapping initial configurations, for which the pseudo-streaming operators 
generate trajectories as well. These trajectories are not time reversal invariant. 

From the explicit form of the pseudo-Liouville operator C, as briefly reviewed in Appendix Al, it is clear that C is 
not odd in the velocity variables, i.e. C c ^ — C, and consequently the trajectory, generated by the streaming operator 
is not invariant under the time reversal transformation. A typical example illustrating the lack of time reversal 
invariance of the hard sphere trajectories, generated by the streaming operators, is a trajectory in T— space starting 
from an unphysical overlapping initial configuration. Here the overlapping pairs separate within a finite time. This 
forward trajectory is not retraced when all velocities are reversed, because the pairs, when meeting again, experience 
a standard hard sphere collision, in which the backward trajectory is deflected from the forward one. 

In the original literature |3Cj it has also been shown that the hard sphere generators obey for t > the commutation 
relation, 

Po e tc = e tZ Po , (II. 7) 

which implies that the time correlation functions are stationary, i.e. depend only a the time difference. In l)A.5)) of 
Appendix A the derivation of (|II.7)I is briefly recalled. In terms of time correlation functions for dynamic variables A 
and B the stationarity relation (|II.2)I also applies, i.e. 

(Ae tc B) Q = J dTAe tZ p B = (Be tc ' A^ , (II.8) 

where the transpose of C is £ e on account of i|II.6|) . Finally we observe that these relations also include conservative 
systems as a special case, in which £ = L on account of (|II.6p and (f II . 3f) . Consequently Eci. (|II.5jl reduces to the 
commutation relation, p$C ■ ■ ■ = Cpo ■ ■ ■ , implying that po is a stationary solution of the Liouville equation, dtp = 
Cp = —Cp. However, Cpo ^ for the Fokker-Planck equation, where Cpo = in the stationary case. 

In summary, the streaming operators for Langevin and hard sphere fluids both lack time reversal invariance, 
and both evolve towards a thermal equilibrium state, and their time correlation functions are stationary. These 
fundamental properties guarantee both the commutation relation l|II.5)l or 1)11.7)) . as well as an identical structure of 
the Green-Kubo formulas for such vastly different systems as hard sphere fluids and Langevin fluids. 

III. ZWANZIG-MORI PROJECTION OPERATOR METHOD 

The goal of this section is to derive for the fluid models, described above, the hydrodynamic equations and Navier- 
Stokes transport coefficients by means of the Zwanzig-Mori projection operator method 0,0. This will be done 
by deriving macroscopic equations for the averages of conserved microscopic densities, or their linear combinations, 
collectively denoted by a*(r), 



a*(r) = {n(r), 5a (r),e(r)}, 



(III.1) 



or by a subset of these densities. Here n(r) = J^. 5(r— rj) is the number density, g(r) = J^. mvi<5(r— r^) the momentum 
density, and e(r) = ^ i ei5(r — r^) the energy density, where the single particle energy ti = \mv 2 + ^ YljM) V{ r ij) or 
the energy of an internal state. In fact this program has been carried out in Ref.[2j for classical fluids with conservative 
interactions. We mostly follow the notation and definitions used in that article. Here we apply the Zwanzig-Mori 
method to the generators of all systems, discussed in the previous sections. 
The general equation of motion for the Fourier-Laplace transform of a l (r,t) is, 

(z - £)a kz = a k , (III.2) 

and Fourier modes are defined as 



flk = / dre-*- r (a(r) - (a(r)) ) 

r DO 

a kz = / dte- zt a k (t), (III.3) 
Jo 

where is a linear combination of conserved densities, which is for small k orthogonal with respect to the inner 
product, 

(ai\a{) = V-VkV k ) = *y + 0(k 2 ), (III.4) 
and the asterisk denotes complex conjugation. This enables us to define a projection operator V = 1 — Vj_, as 

vh * = E = E 4<4i&k>. (ni-5) 



Here (a q |6 k ) = 5 q , k (a^|fe k ) because of translational invariance. Applications of the Zwanzig-Mori method to (|III.2|) 
yields for the projected part ■pa kz , 

[z - V£V - V£P±{z - V±£P±r 1 V±£r}Va^ z = <, (HI.6) 
or in terms of hydrodynamic propagators, Gy(k, z) — (a k |a^ z ), we obtain, 

(111.7) 

3 

The hydrodynamic matrix, Hy = ikftij + k 2 Uij, consists of an equilibrium average, fi^, which is the Euler matrix, and 
a time correlation function, [Ty, which is the Navier-Stokes transport matrix. Its real part, tflHij (k, z), is of 0(k 2 ), 
and represents the dissipative part, and contains the Navier-Stokes transport coefficients of shear viscosity 77, bulk 
viscosity £ and heat conductivity A. For conservative forces ?sHij = kVtij and #lHij — k 2 Uij. For non-conservative 
forces Qij contains real and imaginary parts, as we shall see. By diagonalizing the hydrodynamic matrix for small k 
to 0(k 2 ) one finds the hydrodynamic eigenvalues in the form, z\ ~ ~ikc\ — k 2 D\, and the hydrodynamic modes, a^, 
as special linear combinations of conserved densities. They will be made explicit below. Here the imaginary part of z\ 
contains the propagation speed c\ of that mode (which may be vanishing), and the real part contains the diffusivity 
or dampings coefficient D\ of the mode. We cite the results of equations (44) and (45), obtained in Ref. 0, 

D r , = n/p = Cim ^ {ikn vv {k) + k 2 U m (k,z)} jk 2 

= -£im3i + (C e a2\Vx&V±Cc%)} /& 

Di = -CimuUal\£al) + (£ £ al\V±g z V±£al)} /k 2 

D T = X/(nC p ) = -£imuUal\£al) + (£ e al\Vi_g z V±£al)\ /k 2 . (IH.8) 

Here the symbols D v , Di and Dt represent respectively the transverse, longitudinal and thermal diffusivitics, and 
p = mn is the mass density and C p is the specific heat per particle at constant pressure. Furthermore, £im represents 
the double limit, lirn^o lim/^o ■ I n the equation above g z is the projected resolvent operator, 



g z = v!_(z-v^£VA_Y x Vu 



(III.9) 



which reduces in the small-A:— limit to the standard resolvent, V±(z — C) 1 V±, when acting on projected dynamic 
variables Q. The transformation of the matrix element £/y(k, z) from the representation (a\ s \CP± ■ • • ) in 1111.6(1 to 
the representation (£ e a l k \V± • • ■ ) in l|III.8(l follows directly from the commutation relation pi. 5(1 . valid for Langevin- 
and hard sphere fluids. We also note that the relation applies to systems with dynamics that is time reversal invariant, 
i.e. where C e = —£. 

In summary, the equations as listed in (|III.8|) are valid for all models discussed in this paper. Through this 
formulation we have extended the linear response theory for conservative systems in Ref.jj] to more complex fluids 
like Langevin and hard sphere fluids. 

In order to apply these formulas to different systems the definitions and some properties of the quantities appearing 
in (|III.8|) will be required. Due to the isotropy of simple fluids in thermal equilibrium the (r/i, rjj)— matrix elements 
between shear modes {a^ } are diagonal and independent of the transverse directions (i = 1, 2, ■ • • , d — 1). We further 
recall that the viscosity in isotropic simple fluids is a fourth rank isotropic tensor, determined by two independent 
scalars, i.e. 

Vuf3~/S = V^ajSfjS + SasSfr - (2 / 'd)S a pS lS } + C<W&y<5- (III. 10) 

They are the shear- and bulk viscosity 77 and C respectively, or alternatively the transverse and longitudinal diffusivity, 
D v and Di, determined by, 

Tjxyxy P^t) V 

Vxxxx = pD l = 2r l {l-d- 1 )+C. (III.ll) 

For fluids with conserved densities n(r), <? Q (r) and e(r) — i.e. classical fluids with conservative interactions, and hard 
sphere fluids with impulsive interactions — , the explicit expressions for the hydrodynamic modes are 0, 1^3 1 

< = \/-gk-ki (i=l,2,...,d-l) 




. }j> y -Pk + <ygi \ {a = ±1) 



{e k -/m k }, (111.12) 



nTCp 



where c is the adiabatic speed of sound with c 2 = (dp/dp) s , and pk is the thermodynamic pressure fluctuation, 

*■ (mis) 



dn J \de 

Furthermore {k, k^} (i = 1, 2, . . . , d — 1) is a set of d ortho-normal unit vectors, where the longitudinal direction k is 
taken parallel to the x-axis for convenience. Moreover, g 1 ^ = k • gk is the longitudinal component of gk. Here are 
the shear or transverse velocity modes, a£ (a — ±) the sound modes, and the heat mode. 

It remains to work out the expressions (|III.8fl for the different fluid models, using the appropriate hydrodynamic 
modes in l|III.12|) . in combination with the explicit forms in the pseudo-Liouville operator for these models. 

Equations 1(III.8J) are the main new results of the general part of this paper. We note the latent presence of the 
structure of (|I.2J) . where the instantaneous transport coefficient /ioo oc 5ft (ak|£ak) , and the time correlation functions 
with the different currents, J £ and J, are contained in the matrix J7(k, z). Inserting the currents, defined for small k 
through, 

£a k = -ikJ+0{k 2 ) 

C e at = -iU t + 0{k 2 ), (111.14) 

into pil.8|l reduces this equation to the generic structure (|I.2|I , and we recall that the inner product in (|III.8J| involves 
complex conjugation. Consider first the case of conservative forces. There the dynamics is time reversal invariant, 
C e = —C, i.e. C is odd in the velocity variables, and the contributions of the real part of the Q— matrix necessarily 
vanishes, i.e. /i^ = 0. Moreover, both currents, as defined in (|III.14jl . are equal, J € = J. 

In case the forces are not conservative, i.e. C e ^ — C, then oc 5ftfc~ 2 (a k \Ca k ) is necessarily non-vanishing, 
and approaches a finite value in the small fc— limit. Moreover, it follows from general considerations that J e and J 



are unequal, both for impulsive forces (hard sphere fluids), and for dissipative forces (Langevin fluids). The explicit 
expressions will be discussed in |f|. 

In hindsight the instantaneous transport coefficients /loo can be identified independently of the present analysis. 
In a kinetic theory description, fi^ comes from the irreversible part of the average of the microscopic or mesoscopic 
V-particle current (J)i oc in a state of local equilibrium, i.e. the part of (J)i oc that contributes to the irreversible 
entropy production. In the context of complex fluids this is the part that contributes to the Navier-Stokes transport 
matrix in the hydrodynamic equations. For instance, in the momentum current or stress tensor one finds both in the 
hard sphere fluid |33ll34l|. as well as in Langevin fluids [Tjjj, (Jxy) loc = pu x u y — rjooVxUy. Here the first term, pu x u y , 
is the Euler or reversible part of the current, and the second term is the irreversible part, yielding a contribution 77^ 
to the Navier-Stokes transport coefficients. For the hard sphere fluid the statement above will be explicitly verified 
in Sect. IV. For the Langevin fluids and solids it will be verified in Ref. As the concept of local equilibrium state 
is intuitively clear, one can in principle also estimate fioo for more complex fluids. 

Perhaps the most intuitive understanding of /loc has been given by Hoogerbrugge and Koelman when first 
introducing the dissipative particle dynamics (DPD)-model. The estimate of the shear viscosity of the DPD-fluid 
r\ ~ 7700 by making a continuum (macroscopic) approximation of the Langevin equation of motion. We finally note 
that the instantaneous transport coefficient is also one of the main contributions to the transport coefficients of hard 
sphere- and Langevin fluids at liquid densities. 

The appearance of two different currents J and J c in the Green-Kubo formulas for the viscosity and friction 
coefficient in hard sphere fluids have been derived before in Refs. For an gas of inelastic hard spheres similar 

results have been obtained in Ref. j2||. For a general class of Fokker-Planck equations, satisfying detailed balance, 
local macroscopic evolution equations have been derived in Ref.[2j|- In doing so the expressions, derived for the 
transport coefficients, are also given in the form of Green-Kubo formulas with two different currents J and J e . 
However, the relation between J e and J in that case is more complicated than in the cases discussed in the present 
paper. 

This ends the general part of the paper. The remainder deals with the most important application of the present 
theory, the hard sphere fluid. 



In the application part of this paper we first list for later reference the Green-Kubo formulas in conservative classical 
fluids, and next we work out the details of the Green-Kubo formulas for hard sphere fluids. For classical fluids with 
conservative interactions the explicit forms have been derived in |3j . Here the Liouville operator is odd in the velocities 
because of time reversal invariance. Consequently the diagonal elements n nv and Qtt of the Euler matrix in (|III.8|) 
are vanishing. Computation of the elements of the Navier-Stokes matrix yields for the shear viscosity r\, longitudinal 
viscosity Di and heat conductivity A Q, 



IV. CLASSICAL AND HARD SPHERE FLUIDS 



A. 



Classical conservative fluids 




(IV.l) 



with the subtracted microscopic fluxes, 




Q: 



X 



V ± J X = J, 



X 



p. 



X 




(IV.2) 



Here J xa is the microscopic TV-particle momentum flux or stress tensor, where Fij l<x = —dV(rij)/drij. a with a = 
x,y,...,d is the conservative inter-particle force and V(rij) the spherically symmetric pair potential, {J xx ) = pV 
is the hard sphere pressure in equilibrium, h = (e + p)/n is the enthalpy per particle, and P a = '^ /i mvi a the total 
momentum. Moreover, J a is the total energy flux. The quantity P a h/m, which is the component of J a parallel to 
the conserved quantity P a , is subtracted out in Q x 

B. Hard spheres: Instantaneous transport coefficients 

This section deals with the Cl matrix, (ajJ/Ca^), in which the first major difference in the Green-Kubo formulas 
shows up between systems with time reversal invariance, i.e. C e — —JC, and those without, i.e. C e ^ —L. So, in 
hard sphere systems the Liouville operators necessarily contain even parts in the velocities — apart from possible odd 
parts. The matrix elements of C always have real parts of 0(k 2 ), that contribute to the irreversible or Navier-Stokes 
part of the hydrodynamic equations. These terms are referred to as the instantaneous or high frequency parts of the 
transport coefficients. The appearance of such terms is a very robust feature of all systems with streaming operator 
that lack time reversal invariance. We also recall that the streaming operator exp(tC) in (jiff ,8|l have to be defined 
for all points in phase space, including the physically inaccessible overlapping ones. 

For hard sphere fluids the instantaneous contributions, which are averages in thermal equilibrium, are calculated 
in Appendix A2, and yield results, exact for a d-dimensional hard sphere fluids at general densities, 

rjoa = wdj (d + 2) 

Aoo = wdkB I (2m) = Cyw 

pD ioo = 3zud/(d + 2) 

Coo = vj. (IV.3) 

For completeness we have also listed the instantaneous contribution to the longitudinal diffusivity and the bulk 
viscosity, which have been obtained similarly with the help of (fill. lip . In the equations above we have introduced 

vj = po~ 2 /d 2 tE 

t E = ay / TTf3m/2dbn X , (TV A) 

where the definition of w has been chosen such that it coincides for d = 3 with the notation used in Ref. |34|. and 
tE is the Enskog mean free time, \ = go(r = a + ) is the value of the radial distribution at contact, b = (l/2d)£ldO~ d is 
the excluded volume which equals half the volume of an action sphere with radius a, and fid — 2ir d / 2 /T(d/2) is the 
surface area of a (i-dimensional unit sphere. Alternatively these parameters may be expressed in terms of the equation 
of state for hard spheres through j3p = n(\ + bnx)- 

C. Hard spheres: time correlation functions 

The Navier-Stokes matrix limk-*o U(\l, z) — U(z) gives in the long wave length limit the frequency dependent (z = 
iu>) transport coefficients (memory effects). In the previous subsection we have been dealing with the contributions 
^ocCoojAoo at infinite frequency. We introduce 5C a (t) to denote the inverse Laplace transform of the Navier-Stokes 
part of the hydrodynamic matrix, pU aa (z) = 8C a (z) with a — {i],l,T}. Using the explicit forms l|III.12|l for the 
hydrodynamic modes we obtain from (|III.8Jl the Laplace transforms of the time correlation functions, 

pU m (z) = SC v (z) 

pU u (z) = 8Ci(z) 

uCpUtt(z) = SCt(z) 

The asterisk denotes complex conjugation. From here on only the standard notation of hard sphere kinetic theory 
0, l30l I3TL l32| will be used with £ ~ L + and C e = — L_. To compute the explicit form of the currents we use the 
relation, 



= -|im^ 7 ((/: e 5kx)^^^^k,}o 

= -km^<(£ e 4)P ± 47>±/:e k )o. (TV.5) 



k& k = k5b k = k(b k - S K0 (&o) ) = kSb + 0(k 2 ), 



(IV.6) 



where the fluctuation 5bk = &k — (&k) - Here we have used the relation (&k} = ^k.o (bo) , which is a consequence 
of translational invariance. The replacement of &k by <5&k in l|IV.6(l is convenient for taking the limit k — > 0. It is 
further convenient to choose k parallel to the x-axis. Consider first the longitudinal (a = x) and transversal (a = y) 
momentum current. Then we obtain with the help of l|IV.6|) for small k, 

£-9ka = L + g ka = -ikSJ +xa = -ik(SJ xa + 5J+ xa ) 

C*gi a = = -ik5J_ xa = -tk(SJ k a + SJ v _ xa ). (IV.7) 

By inserting the expressions for the hard sphere generators in IjA.ljl . and using the properties (|A.2|l - l|A.6(l it is 
straightforward to obtain the kinetic part (k) and collisional transfer part (v) of the momentum currents, 

i 

J± xa = ± ^ T±(ij)m(v la r ix + v ja r jx ) 

i<j 

= 2_, mijd \ da(vij ■ d-) 2 a x a a 5(rij — acr). (IV. 8) 

i<j 

The averages of the fluxes are only non-vanishing for a = x, i.e. 

(JL)o = Vn/P, and (J| XX > Q = Vn 2 b x /(3 
(J ±XX ) Q = Vn(l + bn X )/f3 = Vp (IV.9) 

in complete agreement with the result (|IV.2p for smooth potentials. In the transition from the second to the third 
line in l)IV.8(l one may use (|A.8J| and (IA.2|) and express the variables {v^, r^, vj, r 3 } in center of mass coordinates 
{Rij, Gij}, and relative ones {rjj,Vy}, and apply (|A.6|I for conserved variables. Consider first the transverse stress 
fluctuation with a = y. Here we note that V±J± xy = J± xy in l|IV.5|l . as the (xy)-component of the microscopic 
stress tensor has no subtracted parts. The momentum fluxes, J\ xy and JH xy , differ only from one another in the 
constraints on the integral over the d-dimensional solid angle da, denoted by the superscript (— ) or (+). Here (=p) 
requires (T) v ij • <r > 0, and implies that the <r-integration is restricted to the pre-collision hemisphere (— ), or to the 
post-collision hemisphere (+). 

The time correlation function for the stresses for t > can now be written as 

SC v (t) = £ (j- xy e tL + J +xy ) = £ (J +xy e~ tL - J- xy ) = ]T^f(i), (IV.10) 

a, b 

where a, b = {k,v}. The second equality has been obtained by renaming the dummy variables Vj in — Vj and using 
(|A.1|I . Also note that the cross-correlations are equal, SC kv (t) = SC^ k (t), as follows directly from the commutation 
relation l|II.8(l for hard sphere fluids, where C — L + and C e = —L_, i.e. 

{Jx y e + J+xy) = {J+x y e Jxy)o ~ {J-x y e + Jx y )v (IV. 11) 

The last equality has again been obtained by renaming the dummy integration variables to — Vj. The final Green- 
Kubo formula for the shear viscosity in a hard sphere fluid follows then from l|III.8|) . by taking the small- z limit of 
the frequency dependent SC v (z) and combining it with the instantaneous viscosity in (|IV.3j) to yield the equivalent 
representations of the generic structure ljl.2|) . 

/•OO 

V = Voo + {fi/V) / dt { J- xy e tL + J +xy ) . (IV.12) 
Jo 

The hard sphere result for the shear viscosity as well as ljiV.7() - (HV.10(l have been derived already in Ref.0, by 
re-interpreting the Einstein- Helfand formulas in terms of pseudo-Liouville propagators. 

For the heat conductivity the analysis proceeds along similar lines, and the energy flux is calculated from the Fourier 
mode of the energy density e^, in (|III.1|) with ej = \mv 2 , i.e. 



£ek = —ik,J +x = —ik(J x + J+ x ), 



(IV. 13) 



and a relation for J_ x similar to the second line in i|IV.7|) . The subtracted flux is the same as in <|IV.2|) for conservative 
fluids, where the enthalpy per particle for hard spheres is h = h + h v with e = ^dn/(3 and p = (n//3)(l + bnx), and 
where we have used the relation {P X J± X ) = Vn 2 bx/ /3 2 = h v nV//3. The kinetic and collisional parts are then given by 

i 

Q±x = VxJl, = ±imE 1<3 T ± (ij)(v?r ix + v 2 r jx ) - {h-/m)P x 
r(T) 

where Gy = \{yi + Vj-). The time correlation function of the subtracted heat flux is then, 

SCt® = ^ (Q- x e tL +Q +x ) Q = £ SC?(t), (IV.15) 

and the cross-correlations (Aw) and (vk) are again equal (see HIV. 11)1 ). The total heat conductivity in a hard sphere 
fluid is then similarly given by the Green-Kubo formula, 

poo 

A = Aoo + (J3/TV) / dt(Q_ x e tL +Q +x ) Q , (IV.16) 



which has again the generic structure (|f .2(1 . Obviously, this formula can also be derived from the Einstein-Hclfand 
formulas, following the method of Q. 

For the longitudinal diffusivity pDi and the total bulk viscosity C we obtain similarly with the help of (|iV.6J) - (|1V.9J| . 



pDi = P D loo + (P/V) / dt(^_ ra e tL +<5S +;c;c ) 
Jo 

C = Coo + (f3/V) r dt (6S-e tL +6S+) , (IV.17) 



where the subtracted fluxes are, 



SS 



± 



a. 

ss± xx = V±SJ± 2 



{¥) SH-(^) -VV, (TV.IN) 



de J \ dn 

in complete agreement with the subtracted fluxes in l|IV.2(l . Before concluding this section we note that the instan- 
taneous or high frequency transport coefficients are also present in the Enskog theory for dense hard sphere fluids, as 
the local equilibrium averages of the collisional transfer parts in ljIV.8|) and l|IV.14|) for the corresponding currents, i.e. 

{j+xy) loc — —Voc^xUy 
(JlXo = - A ~ V * T 

Y,( J l aa )ioc = P-CooV-u, (IV.19) 

a 

where the bulk viscosity appears as a correction to the local equilibrium hard sphere pressure p, and u is the local 
flow field. As the collisional transfer terms are sums over pairs of particles, the averages contain the pair distribution 
function in local equilibrium, and involves integrations over {vi, v 2 , ri 2 }. They have been calculated in Refs. |33ll34l.l35| 
for the special case of d = 3 as 77 = (3/5)177, A = Cyw and £ = m, and can for instance be found in Eqs.(f 6.f 4,8), 
(16.42,3) and (16. 4f .7)) of Ref.pj] in complete agreement with the results derived here. The explicit expression for 
vj in (|1V.4(I is proportional to n 2 . Consequently at liquid densities the instantaneous transport coefficients are major 
contributions, and account for about half the value of the transport coefficient. 

For comparison we quote the complete prediction of the Enskog theory for the viscosities of a d— dimensional hard 
sphere fluid. For the bulk viscosity one has C = Coo = ro, and for the shear viscosity, 

1 A 2 , \ 2 d 



V=-(l + — 2 bn X ) Vo + j-^. (IV.20) 



Here 770 is the low density Boltzmann value of the shear viscosity |3J| . The expression above shows that the Enskog 
theory gives for the shear viscosity other contributions of order n 2 at high densities, and that the Enskog prediction 
for the instantaneous viscosity 7700 = dw j [d + 2) is exact for all densities. 

V. CONCLUSIONS 

The main result of this paper is a generalization of the standard Green-Kubo or time correlation expressions in 
l|III.8ll for classical fluids with conservative interactions to hard sphere fluids with impulsive interactions, and to 
mesoscopic Langevin fluids with dissipative and stochastic interactions. 

In the classes of fluids above the time dependence of dynamical variables, A(t) = exp[tC}A(0), can be represented by 
streaming operators or generators exp[/j£] when used inside ensemble averages. In conservative systems the dynamics 
generated by the (pseudo-)streaming operators is time reversal invariant, implying Cf = —C (sec Ijll.St ). This property 
is necessary to obtain Green-Kubo formulas of the standard form ljl.l|) . In the remaining classes the dynamics 
generated by the (pseudo-)streaming operator is not time reversal invariant, i.e. the backward streaming operator 
exp(rj£ e ) ^ exp(-/j£) (with t > 0). 

The common criteria, obeyed by the models above are that the systems approach a thermal equilibrium state, 
Pq ~ exp[— f3H], and that the corresponding correlation functions (A(t)B(t')} are stationary, i.e. only depend on 
time differences. The criteria imply that the infinitesimal generators for hard sphere and Langevin fluids obey the 
commutation relation l|II.5|) . i.e. poC — Cpo- This includes in fact conservative systems as the special case with C = C. 

There are two conspicuous differences between Green-Kubo formulas in these three classes of fluids. In the classes 
without time reversal invariant dynamics, i.e. £ c 7^ —£, there is necessarily an instantaneous non-vanishing contribu- 
tion /Ltoo, and in the classes with, i.e. C e = — C, the generator is odd in the velocity variables, and fioo automatically 
vanishes (compare (|IV. L2|) and 1|IV.16|) with (|IV.1J0 . The next difference concerns the microscopic fluxes. In the 
classes with time reversal invariance the Green-Kubo formula in (|I.1|I is a time integral over an a?rfo-correlation func- 
tion ( J(0) J(t)) , containing the same currents. In the cases without, the two currents in (J e e J) Q are unequal. Both 
types of differences apply to the hard sphere fluids, as well as the Langevin fluids. In the former case our results 
confirm and extend the results in |2j. 

A summary of the Green-Kubo formulas for transport coefficients in Langevin fluids can be found in Ref. 5], and 
a detailed analysis will be given in Ref.|jJ, where the results (1111.8(1 will be applied to a Langevin-type solid showing 
heat conduction, and a ditto fluid showing viscosity. 

Finally, to understand the origin of the instantaneous contributions in the case of a hard sphere fluid it is instructive 
to write the shear viscosity (|IV.I2|) in the standard form ijlV.I I) for conservative interactions, where the time correlation 
function for t > is defined as, 

C' v (t) = 2 Voo S(t) + (/3/V) (J^ xy e tL +J +xy ) Q , (V.l) 

and S(t) is a Dirac delta function. It is then of interest to compare this expression with the stress-stress correlation 
function C v (t) in (|iV.IJ) . as recently calculated in Ref. for a conservative fluid with a repulsive power law potential, 
V(r) ~ r~ u (soft spheres) at short times at large values of v. The short time dynamics has been calculated exactly 
for the (very short) mean traversal time r„ ~ l/v, required for particles to traverse the steep part of the repulsive 
potential. It was shown that for asymptotically large v, the short time behavior of C n (t) crosses over to a contribution 
~ S(t), and agrees with the representation l|V.l|) . This also demonstrates that smooth conservative potentials r~ v in 
the singular hard sphere limit as v — > 00, generate a <5(/j)-type short time correlation, in agreement with the present 
results for hard spheres, which have been derived from the very starting point for impulsive interactions. 
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APPENDIX A: HARD SPHERE DYNAMICS 

Al. Pseudo-Liouville and binary collision operators 

The time evolution of any phase space function A(t) — e tc A(0), evolving under the hard sphere dynamics j^lj, can 
be described for positive times (t > 0) by the forward pseudo-streaming operator, exp[ti + ], where C = L + is the 



forward pseudo-Liouville operator for hard spheres. For backward or time-reversed dynamics the pseudo-streaming 
operator is e = e~ tL ~ (t > 0). While true hard sphere dynamics is undefined for physically inaccessible overlapping 
configurations (with any nj < a) , the pseudo-streaming operators are defined for all overlapping and non-overlapping 
configurations. In the latter they generate the true hard sphere dynamics, and in the former they can be defined 
conveniently. Inside averages the pseudo-streaming operators are always preceded by the hard sphere equilibrium 
distribution po ~ Wn exp[— j3K\, where K is the total (conserved) kinetic energy. The overlap function Wn = 1 
for non-overlapping configurations, and Wn = for overlapping ones. Consequently, inside averages the pseudo- 
streaming operators generate the correct forward and backward hard sphere dynamics. For instance, the forward time 
correlation function (Ae tL+ A) Q with t > is equal to the time-reversed one (Ae~ tL - yl) Q . 

In Sec. II we need, apart from the pseudo-Liouville operator £, also its adjoint £, its conjugate £ e with £ = £ c 
according to l|II.6|) . These operators can be expressed in binary collision operators T±(ij) and T±(ij) as follows, 



i<j 

£ e - -£o + J2 T -W = - L - 

i<j 

£ = -L Q + YT-(ij) 

i<j 

£ = Co + Y,T+(ij), (A.l) 

i<j 

where £q — v * ' d/dri, and the binary collision operators satisfy the relations 

T ± =Tl=T T , (A.2) 

as can be verified from the definitions below. The last equality on the first and second line of (|A.1() are the corre- 
sponding notations used in the hard sphere kinetic theory (e.g. fH l32j| ) . The binary collision operators are defined 
as 



T±(ij) = a 4 - 1 J doivtf • &\S(r tj - cjct)^ - 1) 



T±(ij) = o d - x I d&\vij ■ &\ [<5(r y - <j&)b a - 5{rij + ad-)} . (A.3) 

The superscript (=p) denotes the constraint TSij • <r > 0, restricting the ^--integration respectively to the pre- and 
post-collision hemisphere. The 6 CT -operator is a substitution operator, defined as, 

b a Vj = v- = Vi - (v tJ ■ &)& 
b a v, = V ; = Vj + (vy ■ &)&. (A.4) 

The definitions in (|A.lJl imply that the pseudo-streaming operator does not generate a dynamics, that is time- 
reversal invariant, i.e. £ e ^ —£. On the other hand the arguments above imply for hard sphere systems that 
(AB(t)} = (A(—t)B) for t > 0. Expressing this equality in pseudo-streaming operators yields the first equality 
below, i.e. 

r dTp Ae tc B = J dTp Be tc 'A = J dTAe tc ' p B. (A.5) 

In the last equality the adjoint of £ e operator has been introduced. The first equality is based on the properties 
that poey^>{t£) = po exp(t£+) with t > 0, and exp(t£ e )po — exp(— tL_)po with t > generate respectively the exact 
forward and backward dynamics |3f| . So by comparing the utmost left and right expressions we obtain relation l|II.7(l 
for the hard sphere streaming operators in the main text. 

In the body of the paper we also need the property: If a; = {1, Vj, vf} is a collisional invariant, then 



T±(ij)(ai + aj) = 0. 



(A.6) 



A2. Instantaneous transport coefficients 77^ and A, 



Consider first f2 w , which is defined as the instantaneous shear viscosity r/oo, 

fc 2 77oo = -p(a£\Ca£) 

i<3 

where (|A.lfl . (|III.8(I and ljIII.ll) have been used for C and a?, with chosen parallel to the y-axis. Moreover, Cq 
gives a vanishing contribution to l|A.7|) . Calculation of the T + (ij) contribution for small k requires the help of (|A.1(I . 
<|a3J| . (|A~fi)) and (fTVVfj) to show that 

Voc = --y-(^2(viyrix + Vjyr jx )T + (ij)(viyr ix + v jy r jx )) . (A. 8) 

i<j 

From here on the calculations for ijoo are identical to those of (II. 2)- (II. 5) of Ref. 7], where the Einstein-Helfand 
formulas for hard spheres were taken as a starting point. We briefly sketch the different steps. For the hard sphere 
fluid rjoo can be expressed in terms of the hard sphere two-particle distribution function n 2 ipo(vi)(po(vj)x, where 
X = 9o( r = c+) is the radial distribution function at contact and <fo(v) — (Pm,/n) d / 2 exp[— ^f3rav 2 ] is the Maxwellian 
velocity distribution. 

To carry out the subsequent integrations we change to the center of mass, {Rij, Gij = + v j)}j an d relative 
coordinates, {r^ , Vjj = v.j — v,}. The result is 

Voo = -\PP 2 X J dr((g x r y T + (ij)g x r y )), (A.9) 

where g a is the relative velocity, and ((•■•)) denotes a Maxwellian velocity average over all particles involved. Moreover 
the r-integration can be carried out because the operator T + contains S^(r — aa). The remaining integrals are d- 
dimensional generalizations of the collision integrals as appearing in the Enskog theory for hard sphere fluids (see 
Chapter 16.8 of Ref.j^il]). Performing the <r- and g-integrations yields finally for 77^ the result listed in (|IV.3|) of the 
main text. 

The calculations of the instantaneous heat conductivity Aoo runs completely parallel to that of 7700 . Starting from 
(|III.8|) we obtain 



fc 2 Aoo = —nC p (a^\jCa^) 
3_ 
TV 



4r r ((4 " h<) E T +fe')(e k - /m k )>. (A. 10) 



i<] 



By performing the same steps as for rjoo, one sees that the terms involving the enthalpy hn^ do not contribute, and 
one arrives at the result, 

Aoo = -7j^(^2(einx + e j r :jx )T + {ij){e l r ix + ejr jx )) (A. 11) 

i<j 

with ei = \ravl for hard spheres. This expression is completely analogous to (|A.8J) . and a straightforward evaluation 
leads to the expression l|IV.3|) for A,^. 
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